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ABSTRACT 


Ocean modeling is usually constrained by the lack of observed velocity data for 
the initial condition. The diagnostic initialization is widely used to generate velocity data 
as initial condition for ocean modeling. It integrates the model from known temperature 
(T c ), salinity (S c ) and zero velocity fields and holds (T c , S c ) unchanged. After a period of 
the diagnostic run, the velocity field (V c ) is established, and ( T a S a V c ) fields are treated 
as the initial conditions for the numerical modeling. During the diagnostic initialization 
period, the heat and salt ‘source/sink’ terms are generated at each time step. In this 
Thesis, the Princeton Ocean Model implemented to the South China Sea demonstrated 
extremely strong thermohaline sources and sinks generated by the diagnostic 
initialization. Such extremely strong and spatially non-uniform initial heating/cooling 
(salting/freshening) rates in the ocean model may cause drastic change in thermoh aline 
and velocity fields initially (after the diagnostic run). There is a need to overcome such 
problems or find alternative methods as diagnostic initialization is extensively used. 
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I. INTRODUCTION 


A. NUMERICAL OCEAN MODELING 

Numerical ocean modeling is essential to achieve better understanding and 
prediction of the ocean behavior. It integrates hydrodynamic and thermodynamic 
equations numerically with boundary conditions (lateral and vertical) from initial states 
of temperature (T), salinity (S) and velocity. Besides having a good ocean model, it is 
essential to have reliable data for specifying the initial condition in order to achieve 
accurate prediction. Past observations of the ocean have contributed greatly to our 
knowledge of the T and S fields. Relatively reliable climatological (T c , S c ) datasets can 
be obtained from the National Oceanographic Data Center (NODC) or the Navy’s Global 
Digital Environmental Model (GDEM) as initial T and S fields input for the ocean model. 
However, initial velocity field is usually not available due to insufficient number of 
velocity observations and difficulty to obtain such data. Thus, initialization of velocity 
field becomes an important procedure for the ocean modeling. In order to accurately 
predict the ocean behaviors using numerical ocean modeling, a reliable initialization 
process for the velocity field is essential. 

B. PROBLEMS OF THE DIAGNOSTIC INITIALIZATION 

A widely used model initialization for the velocity field is the diagnostic mode. 
This process integrates the model from climatological (T c , S c ), zero velocity fields and 
holds (T c , S c ) unchanged. After a period of about 30 days of the diagnostic run, a quasi¬ 
steady state is achieved and the velocity field (V c ) is established. (T t , S c , Y c ) fields are 
treated as the initial conditions for the numerical modeling. Since initial condition error 
can drastically affect the model’s predictability [Lorenz, 1969; Chu 1999], there is a need 
to ensure that the initial condition data are correctly determined. Since diagnostic mode 
model initialization is being used extensively, there is a need to examine the reliability of 
such an initialization process and determine whether such an approach leads to any non - 
physical phenomenon. 

Chu and Lan [2003] pointed out the problems of the diagnostic initialization and 
found that such an initialization is artificially adding extremely strong heat/salt sources or 
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sinks in the ocean. Their arguments are outlined as follows. The horizontal momentum 
equations for the numerical ocean models are based on 


av av 1 a av 

_ = -v^7 v - W— - k x /V - Vp+ -(K m — ) + H v , (1) 

at dz p dz dz 


where V= (u,v) and w is the horizontal and vertical velocity components respectively. 
V is the horizontal gradient operator. The temperature and salinity equations are 


dT 

dt 


dT d dT 

T -w — +—(K h —)+ H t , 
dz dz dz 


dS 

dt 


-Vg ys-w^- + ^-(K H ^)+H s ., 
dz dz dz 


(2) 

(3) 


where/' is the Coriolis parameter, p the density, p the pressure, and ( K M , K H ) the vertical 
eddy diffus ivity for turbulent mixing of momentum, temperature, and salinity. The terms 
(H v , H t , H s ) represent horizontal diffusion and the subgrid processes causing the local 
time rate of change in (V, T, S). 

The diagnostic initialization is to integrate (l)-(3) under the following conditions: 
T=T C , S =S C , V= 0, at f = 0, (4) 

with T and S unchanged. This is analogous to the process of adding heat and salt 
source/sink terms ( F T , F s ) in (2) and (3). Thus, equations (2) and (3) become 

dT dT d dT , c 

— - T - w—- + — (K h —— )+ H t +F t , (5) 

dt dz dz dz 


dS 

dt 


-Vg y S - w ^- + ±(K H f.) + H s +F s , 
dz dz dz 


and keeping 




( 6 ) 


(7) 
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at each time step. Combining (5), (6) and (7) leads to 


dT d 


dr 


F t = \£7T +w-—-—(K h 


dS d 


dS 


F. = ViPS+w^-~(K h ^-)-H 


(8) 

(9) 


Thus, the heat and salt ‘source/sink’ terms F r and F s are generated artificially at 
each time step during the diagnostic initialization. There is a need to examine how large 
are these source/sink terms generated from the initialization process to determine whether 
these false ‘sources/sinks’ are physical or bearable for numerical modeling. In this thesis, 
POM is implemented for the SCS to evaluate the magnitude of the source/sink terms 
from the diagnostic initialization. Prior to describing the SCS and the model essential 
features, the measures/criteria of the artificial heat and salt ‘source/sink’ terms F T and F s 
have to be established. 


C. CRITERIA FOR STRENGTH OF SOURCE/SINK 


Chu and Lan (2003) had proposed the criteria for the strength of the artificial 
source and sink generated during diagnostic initialization. Based on the SCS condition, 
the maximum annual variability of T, S is estimated to be about 35 °C and 15 ppt 
respectively. Thus, the maximum rates of absolute change of the T, S data are estimated 
as follows: 


dT 

< 35 ° c = oiV 

dS 

dt 

- yr U1 /day ’ 

dt 


15 ppt 
yr 


0.04 PP'/ . 
/day 


( 10 ) 


These values are used as standard measures for ‘source/sink’. Twenty four times of the 
standard measures (10) represent strong ‘source/sink’ and is given as follows: 


dT 


dt 


0.1 


’C, 


Strong 


'hr ’ 


dS 


dt 


0.04 


PPt/ 


I Strong 


'hr 


an 


Ten times of the strong ‘source/sink’ measures (11) represent extremely strong 
‘source/sink’ and is given by: 
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Extremely 

Strong 



dS 

dt \Extre, nely 
Strong 



(12) 


D. PROPOSED RESEARCH IN THIS THESIS 

Although Chu and Lan [2003] found the problem of the diagnostic initialization 
(generation of spurious heat/salt sources and sinks), they did not analyze the uncertainty 
of the initialized velocity field to the uncertainty of horizontal eddy viscosity and the 
duration of the diagnostic initialization. 

The horizontal viscosity is usually taken from the Smagorinsky (1963) formula, 
A M =CAxAy±\AV+(AV) T \ (13) 

where 

i AV + {AVf = (dutdxY+ (dvtdx+du tdy) 1 +(dvldy) 2 ~^, (14) 

C is the horizontal viscosity parameter. 

This thesis is a follow-on work of Chu and Lan [2003] to examine the following 
issues that have not been done earlier: 

1. How long is the duration of the diagnostic initialization needed to get a 
suitable initial velocity field? 

2. How does the uncertainty of the horizontal viscosity parameter C affect 
the artificial heat/salt sources and sinks generated during the diagnostic initialization? 

3. How does the uncertainty of C affect the initial velocity field (Y) derived 
from the diagnostic initiation process? 

4. How is the uncertainty of the velocity field (V) due to the uncertain 
duration of the diagnostic initialization? 

Different from the result presented by Chu and Lan [2003], which was based on 
the Japan/East Sea, this thesis’s area of study is the South China Sea (SCS). The 
Princeton Ocean Model [POM, Blumberg and Mellor, 1987] was implemented for the 
SCS to investigate the physical outcome of the diagnostic initialization with no surface 
and lateral forcing. The NODC annual mean ( T a S c ) data with l°x 1° resolution [Levitus 
et al., 1994] was used. The SCS model was initialized diagnostically for 90 days with 
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four different values (0.05, 0.1,0.2, 0.3) of the horizontal viscosity parameter C. The 60 th 
day velocity field with C = 0.2 was taken as the reference to investigate the uncertainty of 
the diagnostic initialization. 
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II. ENVIRONMEN T OF THE SOUTH CHINA SEA 


The South China Sea (SCS) is the largest marginal sea in the Western Pacific 
Ocean. It includes within its boundaries large shelf regions and deep basins. The deepest 
water is confined to a bowl-type trench between the Philippines and Vietnam, around 
4300 m deep. It includes the shallow Gulf of Thailand and connections to the East China 
Sea (through the Taiwan Strait), the Pacific Ocean (through the Luzon Strait), the Sulu 
Sea (through the Mindoro Strait), the Java Sea (through the Gasper and Karimata Straits) 
and the Indian Ocean (through the Strait of Malacca). It has a total surface area of about 
3.5 x 10 6 km 2 . South of 5° N, the water depth drops to 100m and less as shown in the 
bathymetry chart. 



Longitude (E) 

Figure 1. Bathymetry (m) of the South China Sea (After Chu et al. 2001) 
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All the connecting straits of SCS are shallow except the Luzon Strait, which has a sill 
depth of 2600m. Consequently the SCS is considered a semi-enclosed water body 
(Huang et al., 1994). The complex topography includes the broad shallows of the Sunda 
Shelf in the south/southwest; the continental shelf of the Asian landmass in the north, 
extending from the Gulf of Tonkin to the Taiwan Strait; a deep, elliptical shaped basin in 
the center, and numerous reef islands and underwater plateaus scattered throughout. The 
shelf that extends from the Gulf of Tonkin to the Taiwan Strait is consistently near 70 m 
deep, and averages 150 km in width; the central deep basin is 1900 km along its major 
axis (northeast-southwest) and approximately 1100 km along its minor axis, and extends 
to over 4000 m deep. The Sunda Shelf is the submerged connection between southeast 
Asia, Malaysia, Sumatra, Java, and Borneo and is 100 m deep in the middle; the center of 
the Gulf of Thailand is about 70 m deep. 

The SCS is subjected to a seasonal monsoon system (Wyrtki, 1961). From April 
to August, the weaker southwesterly summer monsoon winds result in a wind stress of 
over 0.1 N/m 2 (Figure 2a) which drives a northward coastal jet off Vietnam and 
anticyclonic circulation in the SCS (Figure 3a). From November to March, the stronger 
northeasterly winter monsoon winds corresponds to a maximum wind stress of nearly 0.3 
N/m 2 (Figure 2b) causing a southward coastal jet and cyclonic circulation in the SCS 
(Figure 3b). The transitional periods are marked by highly variable winds and surface 
currents. 
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25N 


(a) 


(b) 



Longitude (E) Longitude (E) 


Figure 2a and 2b. Climatological wind stress (a) Jun and (b) Dec (after Chu et al. 2001) 



Figure 3. Observational surface circulation: (a) Jun and (b) Dec (after Wyrtki 1961) 


The observed circulation patterns of the intermediate to upper layers of the SCS are 
primarily forced by the local monsoon systems (Wyrtki, 1961), with contributions from 
the Kuroshio Current via the Luzon Strait (or called Bashi Channel), in the southern half 
of the Luzon Strait. The Kuroshio enters the SCS through the southern side of the 
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channel then executes a tight, anticyclonic turn and exits the SCS near Taiwan. Because 
of the kuroshio excursion near the Bashi Strait, occasionally anti-cyclonic rings detached 
from the Kuroshio front would propagate westward into the SCS. This results in a 
persistent southwestward strong current exists near the Dongsha Islands, shown in Figure 
4. 



Figure 4. Kuroshio excursion near the Bashi Strait (after Tomczak 2003) 

Upwelling and downwelling occur off the coast of central Vietnam and eastern 
Hainan. In the north, the waters are cold and saline. The annual variability of salinity is 
small, due to the inflow and diffusion of high salinity water from the Pacific Ocean 
through the Luzon Strait. In the south the tropical conditions cause the waters to be 
warmer and fresher. During the transitions the central region is alternately subjected to 
high and low salinity inflow as the monsoons reverse, resulting in a region of higher 
horizontal gradient and annual variability. Mixed layer depths vary from 30 to 40 m 
during the summer monsoon, and 70 to 90 m (Figure 5) during the winter monsoon with 
variations due to both winds and currents (Wyrtki, 1961). Salinity near the coast of 
major rivers outflow tends to be lower, such as the coast near to Mekong and Pearl rivers, 
shown in Figure 6. 
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Winter 


Summer 



Figure 5. Mixed layer depths of SCS taken from an AXBT experiment 

In winter, the sea surface temperature is generally above 25°C south of 16°N and ranges 
from 20 to 25°C north of 16°N. During summer, sea surface temperature is about 29°C 
south of 16°N and ranges from 25 to 29°C north of 16°N. (See Figure 6 for the sea 
surface temperature plot.) 
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Winter 


Summer 




Longitude 


Longitude 


• Winter: 20-25°C (> 16°N) 

25-27.5°C (< 16°N) 

• Summer: 

25-29°C (> 16°N) 

29-30°C (< 16°N) 

• Pearl & Mekong 
Rivers 


Climatological Salinity 


Climatological Sea Surface Temperature 


Figure 6. Climatological Salinity and Sea Surface Temperature 
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III. MODEL DESCRIPTION 


A. FEATURES OF PRINCETON OCEAN MODEL (POM) 

The Princeton Ocean Model (Alan Blumberg and George Mellor, 1987) is a time- 
dependent, primitive equation numerical model on a three -dimensional grid that includes 
realistic topography and a free surface. The significant attributes of the model are as 
follows: 


1. Sigma Coordinate Model 

It is a sigma coordinate model in that the vertical coordinate is scaled on the water 
column depth. The sigma coordinate equations are based on the following transformation: 


x, y 


y. 


- h 


t = t 


(15) 


H + r i 

where x,y,z are the conventional cartesian coordinates; D = H + r| where H (x, y) 


is the bottom topography and r\(x, y, t) is the surface elevation. Thus, o ranges 
from o=0 at z = t) to o = -1 at z = —H. The sigma coordinate system is 
necessary in dealing with significant topographical variability such as that 
encountered over continental shelf breaks and slopes. Together with the 
turbulence sub -model, the model produces realistic bottom boundary layers, which 
are important in coastal waters (Mellor, 1985). 


2. Free Surface and Mode Splitting 

The model has a free surface and a split time step. For computational efficiency, 
the mode splitting technique [Blumberg and Mellor, 1987] is applied with a barotropic 
time step of 25 seconds, based on the Courant-Friederichs-Levy (CFL) computational 
stability condition and the external wave speed; and a baroclinic time step of 900 
seconds, based on the CFL condition and the internal wave speed. 
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3. Other Significant Attributes of POM 

POM contains an imbedded second moment turbulence closure sub-model to 
provide vertical mixing coefficients. The horizontal grid uses curvilinear orthogonal 
coordinates. The horizontal time differencing is explicit whereas the vertical differencing 
is implicit. The latter eliminates time constraints for the vertical coordinate and permits 
the use of fine vertical resolution in the surface and bottom boundary layers. 

B. SCS MODEL INPUT 

The SCS model contains 125x162 horizontally fixed grid points with 23 o 
levels. The horizontal spacing is ~ 0.179° by 0.175° in the latitudinal and longitudinal 
direction (approximately 20 km resolution). The model domain is from 3.06°S to 
25.07°N, and from 98.84°E to 121.16°E, which encompasses the SCS and the Gulf of 
Thailand. The bottom topography is obtained from the smoothed Naval Oceanographic 
Office Digital Bathymetry Data Base with 5 minutes resolution. The horizontal 
diffusivities are modeled using the Smagorinsky form with the parameter C chosen to be 
0.05, 0.1, 0.2 and 0.3 for this application. No atmospheric forcing is applied to the 
model. 

Closed lateral boundaries, that is, the modeled ocean bordered by land, were 
defined using a free slip condition for velocity and a zero gradient condition for 
temperature and salinity. No advective or diffusive heat, salt or velocity fluxes occur 
through these boundaries. At open boundaries, the radiative boundary condition is used 
with zero volume transport. 

The bottom stress is assumed to follow a quadratic forcing 

x, = (16) 

where p„ =1025 kg m" 3 is assigned as the density of the seawater. V b is the horizontal 
component of the bottom velocity, and D d is the drag coefficient, which is specified as 
0.0025 (Blumberg and Mellor 1987) in the model. 

C. EXPERIMENT DESIGN 

To analyze the impact of the uncertainty of the C to the initialized velocity field, 
one control run and three sensitivity runs of the POM were implemented on the SCS 
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model. The control run was conducted with C = 0.2 while the three sensitivity runs were 
conducted with C = 0.05, 0.1 and 0.3. 

To assess the duration of the initialization and its impact on the velocity field 
under different C-values, the diagnostic model was integrated for 90 days for the four 
cases of C-values. The 60 th day of the model result was used as reference to compute the 
relative root mean square difference of the velocities between day-60 and day-i (i = 60, 

61, 62,.90) to investigate the sensitivity of the initialized velocity field due to the 

uncertainty of the initialization period. 

The POM diagnostic mode was integrated with all three components of velocity 
(u, v', w) initially set to zero, and with temperature and salinity specified by interpolating 
annual mean data to each model grid point. The heat and salt ‘source/sink’ terms F T and 
F s were obtained at each time step. The horizontal distributions of the heat and salt 
source/sink terms at the surface, subsurface, mid-level and near bottom were derived and 
compared to the measures established in (10), (11) and (12). In addition, the horizontal 
mean IF r l and mean IF$1 on the o level with N grid points are calculated in Chapter IV to 
identify the overall strength of the thermohaline source/sink terms generated in the 
diagnostic initiation process. 

The experiment revealed that 30 days were sufficient for the mean model kinetic 
energy per unit mass to reach quasi-steady state under the imposed conditions as 
illustrated in Figures 7 to 10. The thermohaline source/sink terms ( F r , F s ) generated by 
the diagnostic initialization on day-30, day-45, day-60 and day-90 were used to identify 
their magnitudes and sensitivity to the integration period. 
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Figure 7. Mean model kinetic energy per unit mass reached quasi-steady state when POM 
is implemented for the SCS after about 30 days 



Figure 8. IVfean model kinetic energy per unit mass reached quasi-steady state when 
POM is implemented for the SCS after about 30 days 
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Figure 9. IVban model kinetic energy per unit mass reached quasi-steady state when 
POM is implemented for the SCS after about 30 days 



Figure 10. Mean model kinetic energy per unit mass reached quasi-steady state when 
POM is implemented for the SCS after about 30 days 
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IV. RESULT OF SENSITIVITY STUDY - EXTREMELY STRONG 
SOURCE/SINK TERMS IN DIAGNOSTIC INITIALIZATION 

A. HORIZONTAL DISTRIBUTION OF HEAT SOURCE/SINK TERMS 

The criteria (11) and (12) are used to measure the strength of the non-physical 
(artificially created) heat ‘source/sink’ terms generated during the diagnostic 
initialization. The horizontal distributions of F T (°C hr' 1 ) generated by the diagnostic 
initialization at the four levels (near surface, subsurface, mid -level and near bottom) with 
four different C-values show extremely strong heat sources/sinks. The heat sources/sinks 
have various scales and strengths. They reveal small- to meso-scale patterns from the 
surface to the bottom of SCS and include some large-scale patterns near the surface 
(Figures 11 to 14). Near the surface (-10 m) of Taiwan Strait, there is an extremely 
strong heat source present in a large-scale pattern. At this location, the maximum time 
rate of temperature change is 2.217 °C hr 1 which corresponds to an extremely strong 
source of 2642 Wm’ (=p c p F T ) with C = 0.1. Isolated strong heat sources (meso-scale 

pattern) are also present near the surface of Gasper and Karimata Straits, southeast of 
Vietnam and Gulf of Thailand. These features are similar for the four different C-values 
but in general the magnitude of the heat sources decreases when C increases. 

The extremely strong sink is located at the subsurface off Balabac Channel. At 
this location, the extremely strong heat sink reached -3371 Wm 3 (corresponding to F r = 
-2.828 °C hr 1 ) with C = 0.05. Near the northeast of Hainan island and southwest of 
Vietnam (8°N and 104°E), isolated strong heat sinks are present in small to meso-scale 
pattern near the surface. From the subsurface (-100 m) to the mid-level (-1500 m), 
numerous small-scale patterns of strong/extremely strong sources and sinks are well 
mixed and distributed over some areas, such as west of Palawan and northwest of 
Borneo. Near the bottom (-3000 m to -4000 m), the non-physical heat sources and sinks 
are relatively weak. 

In general, when C-value increases, the non-physical sources/sinks weaken, 
nevertheless, they are still above the extremely strong heat source criterion (12) in some 
locations of SCS (as shown in Figure 14 when C = 0.3). It was observed that larger C 
would lead to smaller non-physical heat sources and sinks. However, an excessively 
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large C-value could cause unrealistically strong diffusion in the ocean model and could 
create other adverse consequences to the model results. 
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Figure 11. Florizontal distribution of F T on day-60 with C = 0.05 
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Figure 12. Florizontal distribution of F T on day-60 with C = 0.1 
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Figure 13. Florizontal distribution of F T on day-60 with C = 0.2 
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Figure 14. Florizontal distribution of F T on day-60 with C = 0.3 
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Figure 15. Horizontal distribution of F T near the bottom of SCS (- 4000m) on day-60 
with different C -values 
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B. HORIZONTAL DISTRIBUTION OF SALT SOURCE/SINK TERMS 

The criteria (11) and (12) are used to measure the strength of the non-physical 
(artificially created) salt ‘source/sink’ terms. The horizontal distributions of F s (ppt hr" 1 ) 
generated by the diagnostic initialization at the four levels (near surface, subsurface, mid - 
level and near bottom) with four different C-values show the existence of strong and 
extremely strong salt sources/sinks. 

The non-physical salt sources/sinks have various scales and strengths. They 
reveal small- to meso-scale patterns from the surface to the bottom of the SCS. Near the 
surface, there are also some large-scale patterns (Figure 16 to 19). Near the surface (-10 
m) of Taiwan Strait, there is an extremely strong salinity source present in a large-scale 
pattern. At this location, the maximum time rate of salinity change is 0.422 ppt hr 1 
which is greater than the extremely strong salinity criterion (12) when C = 0.1. Isolated 
strong salinity source (meso-scale pattern) are also present near the surface of Gasper and 
Karimata Straits, southeast of Vietnam and Gulf of Thailand. These features are similar 
for the four different C-values but in general the magnitude of the salt sources decreases 
when Cincreases. 

The strongest salt sink is located at the subsurface of Taiwan Strait. At this 
location, the strong salt sink reached a value of -0.162 ppt hr 1 with C = 0.2 and 0.3. 
Strong salt sinks in meso-scale pattern also present at the surface of Gasper and Karimata 
Straits. From the subsurface (-100 m) to the mid-level (-1500 m), numerous small-scale 
patterns of strong salt sources and sinks are well mixed and distributed throughout some 
areas, such as west of Palawan and northwest of Borneo. Near the botto m (-3000 m to - 
4000 m), the non-physical salt sources and sinks are relatively weak. 

In general, when C increases, the strength of the strong salt sources/sinks 
decreases in magnitude; nevertheless, they are still above the criterion (11) and (12) in 
some locations of SCS (as shown in Figure 19 when C = 0.3). In general, larger C would 
lead to smaller non-physical salt sources and sinks; however, an excessively large C 
value could cause unrealistically strong diffusion in the ocean model and could create 
other adverse consequences to the model results. 
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Figure 16. Florizontal distribution of F s on day-60 with C = 0.05 
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Figure 17. Horizontal distribution of Fs on day-60 with C- 0.1 
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Figure 19. Florizontal distribution of F s on day-60 with C= 0.3 
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Figure 20. Florizontal distribution of F s near the bottom of SCS (- 4000m) on day-60 with 
different C- values 
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C. HORIZONTAL MEAN OF HEAT ‘SOURCE/SINK’ TERM ( M c (I F r I)) 

The horizontal mean IFjl on the <7 level with N grid points is calculated by 

M a (\F T \) = ±-£\Fi\ (17) 

N j=i 

M a (I F r I) is used to identify the overall strength of the heat source/sink terms generated 
in the diagnostic initiation process. Figures 21 to 24 show the temporal evolution of 
M 0 (l/y I) at four different <7 levels: near surface (cr = -0.0125), subsurface ((7 = - 
0.15), mid-level (<7 = -0.5) and near bottom (<7 = -0.95) at different C-values. 




Figure 21. Temporal evolution of M 0 (l/y I) at four different <7 levels with C = 0.05 
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Figure 22. Temporal evolution of M c (I F T I) at four different c levels with C= 0.1 




Figure 23. Temporal evolution of M 0 (I F r I) at four different o levels with C = 0.2 
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Figure 24. Temporal evolution of M 0 (I F r I) at four different cr levels with C = 0.3 


In general, the strength of the absolute mean source/sink decreases across the 
corresponding level when C increases. However, for different C-values, the maximum 
and minimum values of M Q (I F r I) occurred at different levels. When C = 0.05, the 
surface M a (I F r I) increases with time rapidly in the first 5-10 days to reach a peak value 
of 0.033 °C hr 1 and then dips slightly to 0.032 °C hr" 1 . It then oscillates around quasi¬ 
stationary value (0.032 °C hr" 1 ) on the 20 th day onwards. At the subsurface, M 0 (\F r I) 
behaves similar to that at the surface but reaches a higher maximum value of 0.05 °C hr 1 
and a higher quasi-stationary value of 0.0457 °C hr" 1 . At the mid-level, M 0 (I F r I) 
increases with time rapidly in the first 5-10 days to 0.035 °C hr" 1 and then decreases to 
0.033 °C hr 1 on day-20. It then oscillates around 0.033 °C hr 1 from day-20 onwards. 
Near the bottom, M c (I F T I) behaves the same way but with a higher maximum value of 

0.037 °C hr 1 and dips to 0.034 °C hr 1 on the 20 th day and then oscillates around 0.034 
°C hr' 1 . M 0 (\F t I) generally increases with time rapidly in the first 5-10 days and then 
oscillates around quasi-stationary from day-20 or day-30 onwards for the rest of C-values 
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(C = 0.1, 0.2 and 0.3). For C = 0.2, the surface M c (IF) I) increases to 0.025 °C hr 1 in 
the first 5-10 days, then oscillates with decreasing amplitudes until day-30 and finally 
oscillates around 0.021 °C hr' 1 with small amplitudes. The fluctuation patterns from 
subsurface to near bottom are similar to that at the surface except with different 
maximum and quasi-stationary values. These peak values and quasi-stationary values are 
considered large as they represent the average over the entire c level. 

From the depth profile plots of M c (IF) I) on the 30 th day, 45 th day, 60 th day and 
90 th day shown in Figures 25 to 28, when C = 0.05 (Figure 25), the quasi-stationary 
M a (IF) I) has a maximum value at the subsurface (near cr = -0.25) which corresponds 
to a maximum heating rate at the subsurface. The minimum M a (I F r I) is at the mid- 
level (near c = -0.55). When C= 0.1 (Figure 26), the quasi-stationary also M c (\F r I) 

has a maximum value at the subsurface (near c = -0.25) but the minimum value 
occurred at the surface (near cr = -0.0125). Bigger value at the subsurface indicates a 
greater spurious heat sources and sinks problem at this level. However, when C = 0.2 
and 0.3, maximum M a (IF) I) occurred near to the bottom at cr = -0.85 (Figures 27 and 
28). Bigger value near the bottom indicates a greater heat sources and sinks problem at 
this level. The minimum M 0 (IF) I) occurred near the surface. 

From the sensitivity experiments conducted with four different C-values, the 
smallest quasi-stationary value of M c (I F r I) is 0.019 °C hr 1 which is near the SCS 
surface and it occurred when C = 0.3 (Figure 28). The maximum quasi-stationary value 
of M 0 (IF) I) is 0.05 °C hr 1 and it occurred at the subsurface of SCS when C = 0.05 
(Figure 25). 
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Figure 25. Depth profiles of M c (I F T I) with C = 0.05 
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Figure 26. Depth profiles of M a (I F T I) with C= 0.1 
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Figure 28. Depth profiles of M c (I F T I) with C= 0.3 
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D. HORIZONTAL MEAN OF SALT ‘SOURCE/SINK’ TERM ( M 0 (I F $ I)) 

Based on (9) derived in Chapter I, the horizontal mean F 5 I on the cr level is 
calculated by 

M a (IF s l)=-^fV/l, (19) 

A j= 1 

This is used to identify the overall strength of the salt source/sink terms. Figures 29 to 32 
show the temporal evolution of M a (IF 5 I) at four different cr levels: near surface (0 = - 

0.0125), subsurface (o = -0.15), mid-level (c = -0.5) and near bottom (0 =-0.95) at 
different C-values. M a (\F S I) increases with time rapidly in the first 5-10 days and then 
oscillates around quasi-stationary values for the difference cases of C- values. In general, 
the strength of the absolute mean source/sink decreases across the corresponding level 
when C increases. 

When C = 0.05, M c (I F s I) increases with time rapidly in the first 5-10 days to 
reach a peak value of 0.0137 ppt hr 1 (surface) and then dips slightly to oscillate around 
quasi-stationary values on the 20 th day at around 0.0135 ppt hr" 1 . At the subsurface, 
M a (\F S I) again increases with time rapidly in the first 5-10 days to a lower maximum 

(0.0071 ppt hr' 1 ) and then decreases to 0.0068 ppt hr 1 on day-30. It then oscillates 
around 0.0068 ppt hr 1 from day-30 onwards. At the mid-level and near the bottom, 
M 0 (\F S I) behaves similar to that at the surface and subsurface but with a lower 

maximum value and a lower quasi-stationary value. When C = 0.1, 0.2 and 0.3, 
M a (\F S I) behaves similar to that of C = 0.05 but the maximum value and the quasi- 

stationary values decrease when C increases across the corresponding c> level (See 
Figures 29 to 32). 

From the depth profile plots of M c (I F s I) on the 30 th day, 45 th day, 60 th day and 

90 th day shown in Figures 33 to 36, when C = 0.05 (Figure 33), the quasi-stationary 
M a (I F s I) has a maximum value at the surface (near c = -0.0125) which indicates a 

greater artificial salt source/sink problem near the surface. The quasi-steady value of 
M a (I F s I) decreases as the depth increases. The minimum value of M a (I F s I) occurred 

at the bottom for all the four cases of C-values. The smallest quasi-steady value of 
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M a (I F s I) is 0.005 ppt hr 1 near the SCS bottom and it occurs when C = 0.3 (Figure 36). 
The maximum quasi-steady of M G (I F s I) is 0.0135 ppt hr 1 at the surface of SCS when C 
= 0.05 (Figure 29). Bigger M a (I F s I) at the surface indicates a more serious spurious salt 
sources and sinks problem at this level. 
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Figure 29. Temporal evolution of M c (IF^ I) at four different c levels with C= 0.05 
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Figure 30. Temporal evolution of M a (\F S I) at four different cj levels with C = 0.1 
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Figure 31. Temporal evolution of M c (I F s I) at four different o levels with C = 0.2 
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Figure 32. Temporal evolution of M Q (i h\. I) at four different o levels with C = 0.3 
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Figure 33. Depth profiles of M G il F\. I) with C= 0.05 
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Figure 34. Depth profiles of M G (I F s I) with C= 0.1 
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Figure 35. Depth profiles of M n (I F s I) with C= 0.2 
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Figure 36. Depth profiles of M G (I F s I) with C= 0.3 
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V. UNCERTAINTY OF THE DIAGNOSTICALLY INITIALIZED 
VELOCITY FIELD DUE TO THE UNCERTAIN 
HORIZONTAL VISCOSITY 

A. HORIZONTAL VELOCITIES ON DAY-60 OF THE DIAGNOSTIC 

INITIALIZATION 

The horizontal velocities on the 60 lh day of the diagnostic initialization for the 4 
levels (near surface, subsurface, mid-level and near bottom) for each of the four C-values 
are shown in Figures 37 to 40. In general, the surface and subsurface SCS circulation 
heads southward near the east coast of China and Vietnam in an anti-cyclonic pattern for 
all the four cases of C-values. However, such a velocity field could have a Relative Root 
Mean Square Difference as large as 60%. This large uncertainty on the velocity fields is 
shown in the next section. There is also an anti-cyclonic eddy-like structure centered at 
(14°N, 117°E) which appeared from the surface down to the bottom (- 3000m) of SCS. 
Near the bottom of SCS, this anti-cyclonic eddy-like structure is more pronounced when 
C is small (C = 0.05) and appeared to be in a random pattern when C increased to 0.3. 
The anti-cyclonic circulation pattern is contained within the SCS basin and appeared to 
have no exchange with the open boundaries. This is consistent with the model set-up that 
there is zero volume transport at the open boundaries. 

At the subsurface when C = 0.05 and 0.1, there seems to be a strong jet headed 
towards the northwest coast of Borneo. This phenomenon is not apparent when C = 0.2 
and 0.3. From the horizontal velocity plots, there is no obvious relationship between the 
diagnostically initialized velocity field and the C-value. 
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Figure 37. Horizontal velocities on the 60 th day (C=0.05) 
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Figure 38. Florizontal velocities on the 60 lh day (C= 0.1) 


45 
















































































































locffc iaf* liA. iiA laA infs atfE 

Lanjjtmli fdapthi •UQOm) 


IMft 

L ^ Hpd n e»yttr -3900M) 


Figure 39. Horizontal velocities on the 60 th day (C = 0.2) 
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Figure 40. Horizontal velocities on the 60 lh day (C = 0.3) 









B. RELATIVE ROOT MEAN SQUARE DIFFERENCE OF VELOCITIES 


As highlighted in Chapter IV, one of the objectives is to demonstrate how the 
uncertainty of C-value affects the velocity field (V) derived from the diagnostic initiation 
process. To study the impact caused by the uncertainty of C-values, four different C- 
values (0.05, 0.1, 0.2 and 0.3) were used in the POM to conduct four numerical 
experiments. C = 0.2 was used as a control run and the depth-dependent Relative Root 
Mean Square Difference of the horizontal velocity ( RRMSDV) and of the vertical velocity 
(RRMSDW) between the control run and sensitivity runs (C = 0.05, 0.1 and 0.3) were 
computed, 



(17) 


(18) 


where C = 0.05, 0.1 or 0.03. M< is the number of grid points along the east-west 
direction and M, is the number of grid points along the north-south direction for the 
SCS. k= 1 to 23 is the number of a-level. 

Both RRMSDV (k, Cj and RRMSDW(k , C) increase with time rapidly in the first 5- 
10 days and then oscillates around quasi-stationary values for the difference cases of C- 
values. The largest value is between C = 0.05 and C = 0.2 (control run). The value of 
RRMSDW(k, Cj is much larger than the value of RRMSDV(k, C) for C = 0.05, 0.1 and 
0.3. This could due to the relatively smaller magnitude and larger uncertainty of the 
vertical velocity in the SCS. The vertical profile of RRMSDV(k, Cj has a maximum value 
at the mid-level for the different cases of C on day-30, 45, 60 and 90 of the diagnostic 
run. This feature indicates a strong variation of the horizontal velocity in the mid -level of 
SCS. The vertical profile of RRMSDW(k, Q decreases from the surface to the bottom 
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indicates the vertical velocity’s variation decreases from the surface to the bottom. There 
is a decrease in the rate of decrease of RRMSDW(k, Q from the subsurface to the bottom 
for all the three sensitivity runs (C = 0.05, 0.1 and 0.3) when compared with surface to 
subsurface. Thus, the variation of vertical velocity is expected to be smaller from the 
subsurface to the bottom when compared with the surface to the subsurface. Both 
RRMSDV(k, C) and RRMSDW(k , C ) for different C values and o levels are plotted in 
Figures 41 to 52 for illustration. 

Figure 41 shows that RRMSDV(k, 0.05) increases with time rapidly in the first 5 
days to reach 0.6. It then oscillates between 0.6 and 0.8 from the surface to the mid -level. 
Near the SCS bottom, RRMSDV(k, 0.05) oscillates between 0.5 and 0.6. RRMSDV(k, 
0.05) decreases with depth from mid-level to the bottom. The values of RRMSDW(k, 
0.05) and RRMSDW(k, 0.1) are greater than 1 from the surface to the subsurface. 
RRMSDWik, 0.05) decreases to about 0.7 near the bottom (Figure 42). The uncertainty of 
C-value has a significant effect on the velocity field (V) derived from the diagnostic 
initiation process. 
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Figure 43. Depth profiles of RRMSDV(k, 0.05) at different day 
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Figure 44. Depth profiles of RRMSDW(k , 0.05) at different day 
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Figure 51. Depth profiles of RRMSDV(k, 0.3) at different day 



Figure 52. Depth profiles of RRMSDW(k , 0.3) at different day 
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VI. UNCERTAINTY OF THE DIAGNOSTICALLY INITIALIZED 
VELOCITY FIELD DUE TO THE UNCERTAIN LENGTH OF 
DIAGNOSTIC INTEGRATION 

Another unsolved problem is how long the diagnostic integration is needed to 
obtain a suitable initial velocity field (Vc). Figures 1-4 showed that after a period of ~ 30 
days of diagnostic run, a quasi-steady state is achieved. Thus, the 60 th day has been 
selected as reference field to compute RRMSDV(t) and RRMSDW(t) from day-60 to day- 
90, 



(19) 


( 20 ) 


where t = 60, 61.. ..90, M is the number of o-level. 

Figures 53 and 54 show that RRMSDV (t) and RRMSDW(t) fluctuates irregularly 
with time from day-60 to day-90. The value of RRMSDV(t) increases with time rapidly 
from day-60 to day-65 and then oscillates around quasi-stationary values for all the C- 
values. The maximum quasi-stationary of RRMSDV(t) is above 0.6 when C = 0.05. 
When C increases, the quasi-stationary value of RRMSDV(t) decreases. The maximum 
quasi-stationary of RRMSDV(t) is about 0.09 when C= 0.03 (Figure 53). 

The value of RRMSDW(t) increases with time rapidly from day-60 to day-65 and 
then oscillates around quasi-stationary values. The maximum quasi-stationary of 
RRMSDW(t ) is about 0.5 when C = 0.05. When C increases, the quasi-stationary value of 
RRMSDW(t) decreases. The maximum quasi-stationary of RRMSDW(t) is about 0.08 
when C = 0.03 (Figure 54). These oscillations slow down when C increases. Both the 
RRMSDV(t) and RRMSDW{t ) fluctuate irregularly with time (even though quasi-steady 
state has been reached). The uncertainty of the diagnostic integration period affects the 
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uncertainty in the initialized velocity field Vc significantly and in an unpredictable way. 
Therefore, it is difficult to determine the time period that the initial velocity field (Vc) is 
most appropriate. 



Figure 53. Temporally varying of RRMSDV(t) with four different C-values 
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Figure 54. Temporally varying of RRMSDW(t ) with four different C-values 
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VII. CONCLUSIONS 


Diagnostic mode initialization is widely used to obtain the initial velocity field 
(V c ) required for the prognostic model run. Together with the climatological temperature 
and salinity (T c , S c ), (T c , S c , V c ) fields are treated as the initial conditions for the 
numerical modeling. However, POM implemented for the SCS revealed that extremely 
strong thermohaline source/sink terms were generated for C = 0.05, 0.1, 0.2 and 0.3. 
Such non-physical phenomenon was also observed in the Japan/East Sea when diagnostic 
initialization was conducted using POM by Chu and Lan [2003]. 

In this thesis, we saw that the strength of the mean heat source/sink, M Q (\F r I), 
fluctuates with depth irregularly from the surface to the bottom. The fluctuation also 
varies with C. At different C-values, the maximum and minimum quasi-stationary values 
of M c (i Fj I) occurred at different levels. In general, M 0 (I F T I) increases with time 
rapidly in the first 5-10 days and then oscillates at different quasi-stationary values for the 
different C-values. The strength of the absolute mean source/sink decreases across the 
corresponding level when C increases. The horizontal distributions of heat source/sink 
term, F r (°C hr" 1 ), at the four levels (near surface, subsurface, mid -level and near bottom) 
at four different C-values show extremely strong heat sources/sinks generated by the 
diagnostic initialization. The heat sources/sinks have various scales and strengths. They 
revealed small- to meso-scale patterns from the surface to the bottom of the SCS but 
there are some large-scale patterns near the surface of the SCS. In general, when C 
increases, the strength of the strong and extremely strong sources/sinks decreases in 
magnitude; nevertheless, they are still above the criteria (11) and (12) in some locations 
of the SCS. Near the bottom (-3000 m to -4000 m), the SCS basin is less affected by the 
spurious heat sources and sinks. Although larger C would lead to smaller spurious heat 
sources and sinks, an excessively large C-value could cause unrealistically strong 
diffusion in the ocean model and could create other adverse consequences to the model 
result. 

The strength of the mean salt source/sink, M c (I F s I) decreases with depth from 
the surface to the bottom. It increases with time rapidly in the first 5-10 days and then 
oscillates around quasi-stationary values for the difference cases of C-values after day-20. 
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These features are consistent for the four different C-values but the strength of M a (I F s I) 

decreases across the corresponding (7 level when C increases. The horizontal 
distributions of the salt source/sink term, F s (°C hr" 1 ), at the four c levels (near surface, 
subsurface, mid-level and near bottom) at four different C-values generated strong and 
extremely strong salt sources/sinks from the diagnostic initialization. They revealed 
small- to meso-scale patterns from the surface to the bottom of SCS. In addition, there 
are some large-scale patterns near the surface of the SCS. Greater salting rate occurred at 
the surface indicated a more serious spurious salt sources and sinks problem at this level. 

The effect of the uncertainty of C on the velocity field generated by the diagnostic 
initialization process was analyzed through the use of Relative Mean Square Error. It 
was observed that RRMSDV(k, Q and RRMSDW(k, C) increase with time rapidly in the 
first 5-10 days and then oscillate around quasi-stationary values for the difference cases 
of C-values. The largest uncertainty is between C = 0.05 and C = 0.2 (control run) and 
the smallest is between C = 0.3 and C = 0.2 (control run). The uncertainty of the 
diagnostic integration period affects drastically the uncertainty in the initialized velocity 
field Vc. These were demonstrated by the RRMSDV (t) and RRMSDW(t) which fluctuate 
irregularly with time from day-60 to day-90. Thus, it is difficult to determine the time 
which will yield the most appropriate Vc. 

When the prognostic integration starts, F T and F s are immediately removed from 
(5) and (6) (shown in Chapter I). The extremely strong and spatially non-uniform initial 
heating/cooling (salting/freshening) rates are introduced into the ocean models and cause 
drastic changes in the thermohaline and velocity fields initially (after the diagnostic run). 
This effect is especially significant in the deep layer below the thermocline and halocline. 
Since diagnostic initialization is extensively used, there is an immediate need to develop 
a check-up algorithm to monitor the strength of the ‘source/sink’ terms during the 
diagnostic initialization process; otherwise, alternative methods would have to be used to 
derive Vc. 
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APPENDIX A. HORIZONTAL DISTRIBUTION OF HEAT 
SOURCE/SINK TERMS 


Horizontal distribution of F T on day-30 with C = 0.05 
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Latitude Latitude 


Horizontal distribution of F T on day-90 with C = 0.05 
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Horizontal distribution of F T on day-30 with C = 0.1 
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Latitude Latitude 


Horizontal distribution of F T on day-90 with C= 0.1 
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Latitude latitude 


Horizontal distribution of F T on day-30 with C = 0.2 
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Horizontal distribution of F T on day-90 with C = 0.2 
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Horizontal distribution of F T on day-30 with C = 0.3 
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Horizontal distribution of F T on day-90 with C = 0.3 
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APPENDIX B. HORIZONTAL DISTRIBUTION OF SALT 
SOURCE/SINK TERMS 


Horizontal distribution of F s on day-30 with C = 0.05 
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Horizontal distribution of F s on day-90 with C = 0.05 
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Horizontal distribution of F s on day-30 with C = 0.1 
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Horizontal distribution of F s on day-90 with C = 0.1 
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Latitude Latitude 


Horizontal distribution of F s on day-30 with C = 0.2 
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Latitude Latitude 


Horizontal distribution of F s on day-90 with C = 0.2 
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Latitude Latitude 


Horizontal distribution of F s on day-30 with C = 0.3 
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Latitude Latitude 


Horizontal distribution of F s on day-90 with C = 0.3 
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